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ABSTRACT 

Recent studies have shown that baroclinic vortex amplification is strongly dependent on certain 
factors, namely, the global entropy gradient, the efficiency of thermal diffusion and / or relaxation as 
well as numerical resolution. We conduct a comprehensive study of a broad range and combination 
of various entropy gradients, thermal diffusion and thermal relaxation time-scales via local shearing 
sheet simulations covering the parameter space relevant for protoplanetary disks. We measure the 
Reynolds stresses as a function of our control parameters and see that there is angular momentum 
transport even for entropy gradients as low as ^ — —dlns/dlnr = 1/2. Values we expect in proto- 
planetary disks are between /5 = 0.5 — 2.0 The amplification-rate of the perturbations, F, appears to be 
proportional to and thus proportional to the square of the Brimt-Vaisala frequency(F oc jg^ oc N^). 

The saturation level of Reynolds stresses on the other hand seems to be proportional to This 
highlights the importance of baroclinic effects even for the low entropy gradients expected in proto- 
planetary disks. 

Subject headings: accretion, accretion disks, circiraisteUar matter, hydrod5mamics, instabilities, turbu- 
lence, methods: numerical, solar system: formation, planetary systems 



1. INTRODUCTION 

Angular-momentum transport and turbulence are 
important issues concerning protoplanetary disks. 
Magneto-hydrod5mamic turbulence brought about by 
the magnetorotational instability (MRI, Balbus & Haw- 
ley 1991), is a reliable way to achieve a sufficient 
angular-momentum transport and with this also an ac- 
cretion rate fitting observations (Andrews et al. 2009) 
and playing an important role in planet formation (Jo- 
hansen et al. 2007; Lyra et al. 2008; Dzyurkevich et al. 
2010; Flock et al. 2011; Uribe et al. 2011; Johansen et al. 
2011). However, for MRI to be active the gas has to be 
sufficiently ionized. This is only the case in the outer 
regions, upper layers of the disk, and in regions close 
to the star The other parts of the disk are too cold and 
dust-rich for sufficient ionization and thus the magnetic 
fields cannot couple to the gas. Because of this, the MRI 
carmot operate in this region, which is therefore called 
"dead zone" (Gammie 1996; Turner & Drake 2009). 

Since the precise ionization structure is still under 
debate (Turner & Drake 2009) as is the interplay be- 
tween active and dead-zones (Lyra & Mac Low 2012) 
we want to assess the precise hydrod5mamic behavior 
of dead zones, because accretion has to proceed through 
it somehow and it is where planets form. Therefore it 
is of interest to study purely hydrod5mamic turbulence 
in circumstellar disks. Klahr & Bodenheimer (2003) 
found such a hydrod5mamic instability creating vortices 
in three-dimensional radiation hydrod5mamical simula- 
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tions of baroclinic disks, e.g. with a radial entropy gradi- 
ent and thus vertical shear, which they assumed to be a 
kind of baroclinic instability (BI) modified by the Keple- 
rian shear profile. Observed protoplanetary disks have 
a non-zero radial entropy gradient j6 = —d\ns/d\nr, 
where s is the entropy and r the radial distance to the 
star With ^ — q — (72D — 1) PE/ where q = —dlnT/dlnr 
and pz — — dlnE/dlnr are the temperature surface den- 
sity gradient respectively and 720 is the 2D adiabatic 
index, we see that disks that fulfill < q/ (72D — 1) 
indeed have a negative entropy gradient with values 
from Andrews et al. (2009) of q ^ 0.3 - 0.5 and = 0.9. 
Therefore protoplanetary disks are not barotropic but 
rather baroclinic which means that planes of constant 
pressure and constant density are misaligned, creating 
a thermal wind, e.g. vertical shear. In a linear stability 
analysis that followed (Klahr 2004) it was shown that 
this instability can only be of non-linear nature (see also 
Cabot 1984; Knobloch & Spruit 1986). 

Thermal relaxation turned out to be crucial when Pe- 
tersen et al. (2007a,b) studied baroclinic vortex amplifi- 
cation using an incompressible approximation. In fact 
thermal relaxation or diffusion, besides the entropy gra- 
dient, are key ingredient to establish baroclinic feedback 
that keeps the instability e.g. vortices in baroclinic disks 
growing. 

While both effects e.g. the baroclinic instability an.d 
baroclinic vortex amplification are a result of the su- 
peradiabatic radial stratification of a disk they are not to 
be confused. An operating linear baroclinic instability 
(compare Cabot 1984; Knobloch & Spruit 1986) would 
be able to create vortices in disks from infinitesimal per- 
turbations, whereas the baroclinic vortex amplification 
deals with the growth of existing vortical perturbations, 
for which Lesur & Papaloizou (2010) used the term "sub- 
critical baroclinic instability" (SBI). 

The occurrence of a classical BI in the disk in its geo- 
physical definition is stiU imder debate and shall be dis- 
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cussed elsewhere. There are three possibilities: 1.) there 
is a classical BI working in protoplanetary disks creating 
the initial vortices, 2.) there is an other instability op- 
erating (see the discussion in Klahr 2004) for instance 
creating vortices via Kelvin-Helmholz instability from 
vorticity maxima in sheared waves of baroclinic disks 
or 3.) small vortical perturbations are triggered from 
other effects, e.g. waves from the MHD active region of 
the disks or maybe from the waves emitted by vortices 
at other radii. In any case the vortices are then grow- 
ing as described by the BVA until they reach a sufficient 
size to influence the evolution of the disk, and this is the 
physics being subject of the present paper. 

Recently, Lyra & Klahr (2011) have examined the in- 
terplay of baroclinic vortex amplification and MHD. 
They found that as soon as magnetic fields are coupled 
to the gas, the MRI takes over and thus superseeds vor- 
tices which were previously amplified by vortex ampli- 
fication. This is evidence that the vortex amplification is 
a phenomenon restricted to the dead-zone. 

All the above mentioned (lower resolution) studies 
had to apply entropy gradients 2-4 times stronger than 
to be expected in protoplanetary disks (Andrews et al. 
2009, Klahr 2013 submitted) to drive BVA. We show in 
the current paper, through high resolution runs that re- 
alistic entropy gradients in protoplanetary disks are suf- 
ficient for BVA. 

Recently Paardekooper et al. (2010) have investigated 
the effect of radial vortex migration. They discov- 
ered that vortices migrate quickly radially inward once 
grown to their full size. While this effect will be of ma- 
jor importance to understand the life-cylce of a vortex, it 
plays a weaker role for the small /still growing vortices 
in the present paper. Of course migration will influence 
the effective angular momentum transport generated by 
the vortices via the emission of waves, but this is be- 
yond the scope of 2D local simulations as in our study. 
We shall return to vortex migration and have a better es- 
timate for angular momentum transport once we return 
to global simulations. 

We carry out local, compressible shearing sheet sim- 
ulations at various resolutions. We show that as we go 
to higher resolutions one can excite the nonlinear insta- 
bility and achieve Re5molds stresses with the low en- 
tropy gradients deduced for observed accretion disks. 
We conduct an extensive parameter study for entropy 
gradients (/S), resolution, thermal cooling (Tcqoi) and dif- 
fusion times (x^iff) respectively. Section 2 gives a brief 
overview of the physical background of the instability. 
In Section 3 we present the numerical setup of our sim- 
iilations. In Section 4 we examine the amplification and 
decay-times of values such as enstrophy cv^ = {V x u)l 
and a-stresses. Here a = {pUxUy{qpQ)^^) with p being 
the gas density, u the gas velocity, q — 1.5 the shear pa- 
rameter, and po the initial mean pressure. We also ana- 
lyze the saturation values, e.g. how quantities like the 
entropy gradient, cooling processes in the disk or the 
size of the simulated domain influence the strength of 
angular momentum transport. Finally we summarize 
our results and give a conclusion in Section 5. 
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Vorticity is conserved in quasi-incompressible 
barotropic simulations, but in flows with density 
and pressure as independent quantities vorticity is 
produced via the so called baroclinic term 

^ = Vx (^-Ivp) =^V^)x Vpoc/Sa^p. (1) 

Here p is the gas density, p the gas pressure, and j6 is 
the global radial entropy gradient. The ground state 
of a disk is geostrophic, e.g. all centrifugal forces and 
gravity are in balance with the strictly radial pressure 
gradient. If an entropy perturbation is introduced with- 
out perturbing the pressure, then this entropy pertur- 
bation will efficiently create vorticity in the presence of 
the global entropy and pressure gradients. This effect is 
basically radial buoyancy because of superadiabatic ra- 
dial stratification^. Indeed the radial Brunt-Vaisala fre- 
quency (Tassoiil 2000) 

f^) (2) 

jp or or \p^ J 

is imaginary, which would lead to radial convection. 
However, shear stabilizes non-axisymm.etric modes and 
for the dynamic stability of the axisymmietric system the 
Solberg-Hoiland criterium (Tassoul 2000; Rudiger et al. 
2002) 

1 df 1 

dp / dj^ ds df- ds\ 
dz \dRdz ~ 'dzdRJ ^ 

has to be considered. If one re- writes Eq. (3) for local 
approximation (see e.g. Balbus & Hawley 1998) the sta- 
bilizing action of the specific angular momentum shows 
up as the value of Oort's constant in the Coriolis term. If 
also the vertical stratification in velocity is taken under 
consideration, as it will occur in real three-dimensional 
accretion disks (Fromang et al. 2011), then the combined 
action of radial buoyancy and Coriolis forces lead to a 
thermal wind, e.g. a vertical shear in rotational velocity. 
This is precisely the initial state as baroclinic instability 
in rotating stars and planetary atmospheres. Yet, insta- 
bility in these systems is not obstructed by radial shear, 
whereas in a Keplerian disk radial scales would have to 
be on the order of the vertical pressure scale-heigth (H) 
(Knobloch & Spruit 1986) to be linearly unstable with 
respect to baroclinic instability. 

Before we explain the motion of a gas parcel in a vor- 
tex we want to explain the cooling and heating pro- 
cesses in a disk as they proved to be crucial to maintain 
the baroclinic feedback (Petersen et al. 2007a,b). Dust 
particles absorb photons which heats them up. To cool 
they radiate photons in the infrared. This radiation can 
be absorbed by other particles. This happens on a typ- 
ical length-scale. A convenient parametrization for the 

^ Note that similar situations can be found in subadiabatic config- 
urations. In fact, in any non barotropic disk, an entropy perturbation 
will lead to a vorticity fluctuation. But without the global pressure and 
entropy gradient pointing in the same direction these perturbations 
will quickly decay (shear away) as they are lacking the mechanism of 
vortex ampUflcation. 
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diffusion time in our vortex system is T^iff = a^ /K where 
a is the radius of the vortex and K the diffusion constant. 
The diffusion constant can be approached using a flux 
limited diffusion approach as in Kley et al. (2009). There 

K = Ac4flRT^ {p^)~ where A is the flux limiter, c the 
speed of light, the radiation constant, T and p the gas 
temperature and density, respectively and k the opac- 
ity. Since K is constant and the vortex grows T^iff wiU 
change over time. Thermal relaxation is the other pro- 
cess by which dust can deposit heat into the gas. When 
a dust particle has a certain temperature other than the 
equilibrium temperature it will exchange heat with the 
ambient medium until it reaches the backgroimd tem- 
perature again. TcqoI is the time needed to achieve this. 
This time-scale affects vortices of all sizes equally. 

The baroclinic feedback itself was explained in detail 
by Petersen et al. (2007b). A nice description of the 
mechanism can also be found in Lesur & Papaloizou 
(2010). In a baroclinic flow entropy is a function of pres- 
sure and density, s {p,p)- Pressure on the other hand is 
only a function of radius. The vortex interior transports 
high entropy material from small radii to large radii. Af- 
ter thermalization low entropy material is transported 
to small radii. Since the pressure variations, especially 
from weak vortices, are negligible in comparison to the 
global radial pressure gradient and much smaller than 
the azimuthal entropy gradient, pressure can be seen 
as approximately azimuthally constant (Klahr & Boden- 
heimer 2003; Klahr 2004; Petersen et al. 2007a). To keep 
the pressure constant an azimuthal density gradient is 
established, e.g. outflowing material has a lower den- 
sity as inflowing material. Thus the vortex feels the ef- 
fect of differential buoyancy which established the pos- 
itive baroclinic feedback (Eq. (1)). 

If cooling is too fast (short time-scales) then the fluid 
parcel adapts the background temperature slope too 
quickly. The vortex becomes locally isothermal and no 
entropy transport is possible. Conversely, if cooling is 
too slow (long time-scales) then gas will not be ther- 
malized fast enough. The vortex gas becomes adiabatic 
with constant entropy across the vortex. In both extreme 
cases, isothermal or adiabatic, the azimuthal entropy 
gradient across the vortex vanishes. As shown in Eq. 
(1) the vorticity source ceases to amplify the vortex, or 
at least stabilizes it against losses from numerical viscos- 
ity from radiating vorticity perturbations, e.g. Rossby 
waves. Therefore it is important that thermal cooling 
and diffusion times are in the right regime. 

We model both thermal relaxation and thermal diffu- 
sion separately because, dependent on the vortex size, 
either one or the other dominates thermalization. Al- 
ways the process with the shorter time-scale sets the 
heat exchange between vortex and ambient gas. 

3. NUMERICAL SETUP 

Our simulations were conducted with the PENCIL 
CODE^. We use a two-dimensional, local shearing sheet 
approach. We consider a sheet in the mid-plane that co- 
rotates with the co-rotational radius Rq. This is a 2D 
version of the model used in Lyra & Klahr (2011). To 
include the baroclinic term they define a global entropy 

* See http://www.nordita.org/software/pencil-code/ 



gradient f>. Note that in our approximation the gradi- 
ents for entropy (s) and pressure (p) are the same. There- 
fore we do not distinguish between them in our notation 
and call both j6. However, in real disks both may easily 
differ. 

The total pressure ptot = p + p consist of a local fluctu- 
ation p and a time-independent part that follows a large 
scale radial pressure gradient /3 



p^po{r/Ro)~^ 



(4) 



where r is the cylindrical radius. The full set of lin- 
earized equations used in our simulations is 

'^ + ^u-V)p^-pV-u+Mp) (5) 

+ (m ■ V) M = - ^ Vp - 2no (z X u) 



vt 



3^ . /Spo/l 1 
2 Rq \p po 



^x + fv{u,p) (6) 



+ 



/3po u. 

Ro (7- 



(7) 



Here p is the gas density, u is the deviation of the gas 
velocity from the Keplerian value, T the temperature, 
Cj, the specific heat at constant volume and, K the heat 
conductivity. Tthermal diffusion time-scale is denoted 
by Tcooi- The symbol 



V 

m 



5 , (0) 9 



dy 



(8) 



,(0) 



represents the Keplerian derivative where Uy 

-3/2nox. 

For a more thorough derivation of these equations 
and the linearization of the global pressure gradient we 
refer to Lyra & Klahr (2011) and the appendix therein. 

In order to keep the numerical scheme stable we 
add sixth-order hyperdiffusion /d(|0), hjrperviscosity 
fy{u,p), and hyperconductivity /k(s) (Lyra et al. 2008, 
2009; Oishi & Mac Low 2009). 

The radiation processes in the disk are implemented 
through the first (thermal diffusion as an approximation 
for flux limited diffusion of radiation energy density) 
and second (thermal relaxation to mimic heat exchange 
with the surface of the disk and thermal equilibration 
with the irradiation from the central object) terms on the 
right hand side of the entropy equation. As mentioned 
in the last chapter we keep the diffusion coefficient K, 
which is defined as in (Kley et al. 2009), constant and 
define its value via T^iff — H^/K. So if the vortex has a 
radius of H, the pressure scale-hight of the disk, the dif- 
fusion time Tdiff has the value we quote in e.g. Table 1. 
If the vortex is smaller than H relaxation wiU be much 
faster. 

To clarify that it is indeed the global entropy gradi- 
ent that produces the vorticity we take the curl of the 
Navier-Stokes Eq. (6) and assume an eqiiilibrium state. 
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TABLE 1 
Simulation setup and results 
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Ux = 0, and VP = so that 



Vt p^Rq 



dyp. 



(9) 



Here we see that the negative azimuthal density gradi- 
ent across the vortex is the source for vorticity produc- 
tion proportional to the global entropy gradient. 

Shearing sheet simulations with Zeus'' like finite vol- 
ume codes without explicit viscosity, e.g. the TRAMP 
code, have shown a weak amplification of kinetic en- 
ergy for the pure adiabatic case, i.e. infinite cooling 
time (see Klahr 2013 ApJ submitted). This numerical 
artifact does not occur with simulations performed by 
the Pencil Code. See Appendix A for a ID radial 
test/ comparison simulation. 

Initially we apply a finite perturbation in the density 
so that 

p{x,y)=pQ+p' (10) 

with jOo the constant backgroimd density and p' the ac- 
tual perturbation of the form 



P ^po 



Ce-(-/2af 



X . y 



where C describes the strength of the perturbation. We 
perturb the density in a way that jOrms — 5% for ^ — 
1.0,2.0 (runs A-I) and prms = 10% for = 0.5 (runs J- 
P). To achieve a random perturbation we apply an ar- 
bitrary phase (pij between and 1. The initial state is 
non-vortical. Again, this is the identical initial condi- 
tion as used in Lyra & Klahr (2011) as well as the same 
amplitude, C, for simulations with /3 — 2.0, as was used 
in their simulations. 

Note that with this initial perturbation we do not per- 
turb the pressure but the entropy. Thus it is really only 

http://www.astro.princeton.edu/~jstone/ zeus.html 



the term in Eq. (9) that creates the development of non 
laminar flow structtire. 

All our simulations are done in dimensionless code- 
units. So that _Ro = ^"^0 = 1/ 7 = 1-4, and Cs = 0.1, which 
means that H — 0.1. All time-quantities are given in 

27rnQ ^ which is one local orbit at the co-rotational ra- 
dius Rq- 

The individual setups are given in Table 1. The ther- 
mal cooling times and thermal diffusion times are de- 
rived from standard disk models like in Bell et al. (1997), 
also see Klahr 2013 submitted. 

We explored different resolutions in our simulations, 
namely 288^, 576^ and 1152^. The unusual non power 
of 2 resolution comes from our computational platform 
with 6 core processors. Typically we used up to 24 CPUs 
totaling 144 cores for our largest grids. Still we needed 
about 1200 hours per run. The grid covers ±2H around 
Rq in the radial and [0H,16H] in azimuthal direction. 
This leads to an effective resolution of 72 (2882), ^44 
(576^) and 288 (1152^) grid-points per scale hight in ra- 
dial direction and 18 (288^), 36 (576^) and 72 (1152^) 
grid-points per H in azimuthal direction. It is always 
necessary to compromise between resolution and com- 
putational time. Lower resolution simulations are com- 
putationally less expensive but might not resolve the 
necessary scales. 

4. RESULTS 
4.1. Saturation Values and Convergence 
We show the time-developement of a-stresses in 
Fig. 1. The green line shows the resolution of 288^, 
black of 576^ and red 11522 for ^ = 2.0 (top), ,6 = 1.0 
(middle) and /5 = 0.5 (lower panel). In all simulations 
■^diff = '''cool = 10 local orbits. 

We see that for /S = 1.0 and 0.5 and a resolution of 288^ 
the perturbation decays right away. Higher resolution is 
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200 400 600 800 1000 

time [27i/i2„] 



Fig. 1. — Time evolution of a-stresses for tfie tliree different resolu- 
tions of 288^, 576^ and 1152^ with an entropy gradient of /3 = 2.0 (green 
line), /5 = 1.0 (black line) and fi = 0.5 (red line). For all these models 
Tdiff = Tcool = 10 ■ 27r/ Dq. For all resolutions vortex amplification and 
therefore angular momentum transport can be seen for strong entropy 
gradients (/5 = 2.0). For lower entropy gradients higher resolution is 
needed to see the development of vortices. The dashed lines show the 
saturations values {fi = 2.0 and fi = 1.0) and value at the end of the 
simulation {fi = 0.5) respectively. 

required to increase the Reynolds-number of the system 
and have less dissipation on the smaller scales and thus 
excite the instability again. 

We take a stronger initial perturbation for fi — 0.5 than 
for the higher fi. The perturbation in entropy results in 
a perturbation in vorticity. This perturbation is propor- 
tional to fi. For small fi we have to apply a stronger per- 
turbation to get the same effect on the vorticity. How- 
ever, we expect that if we go to even higher resolution 
it is possible to keep the initial density perturbation at 
|Drms = 5% (Petersen et al. 2007). 

If we compare the saturation values of runs with dif- 
ferent resolution, we see that they differ by only 10 % 
from one another (see Table 1). 

It is important to note that the instability is excited 
and we measure a-values in the converged runs up to 




1500 



time [27l/Qo] 

Fig. 2. — Time evolution of the a-values and enstrophy for fi = 1.0 
and a resolution of 576^ (run C). The red slope marks exponential 
amplification with a amplification-time t = 70 f^. For larger en- 
tropy gradients (smaller entropy gradients) we get faster (slower) 
amplification- times . 

4 X 10~^ for entropy gradients as low as /3 = 0.5. In fact, 
in Section 4.5 we show that there is only a weak depen- 
dence of a on /5 as a o< fi^-^. Fig. 1 shows that the satura- 
tion values of a do not depend strongly on fi, but as we 
will see in the next section the amplification rates do. 

4.2. Amplification- and Decay Rates 

We analyze the amplification timescales of the vor- 
tices, meaning how fast a vortex grows due to the baro- 
clinic feedback. Thus it is independent of the precise 
shape of the initial condition as long as the amplitude 
is large enough for the given Re3molds number to have 
vortex growth. In fact, the initial strong kick needed to 
get the vortex going decays rather quickly as can be seen 
in e.g. Fig. 1. Here, the a-values start out in the order of 
10~^ then drop to around 10^^ as the initial perturba- 
tion decays. As soon as the baroclinic feedback sets in, 
the values rise again. The timespan that follows is the 
one where we measure the amplification time. 

In analyzing the amplification-rates of the instability 
we find that the initial amplification-rate of the a-stress 
(F (a)), as can be seen in Fig. 2 for run C, can be fit- 
ted as exponential amplification a = dQexp {t / t) with 
T ~ 70/5^^. The proportionality to fi^^ is not what one 
would naively expect from a linear convective or buoy- 
ancy driven turbulence. 

For a linear buoyancy driven turbulence one would 
expect an amplification rate proportional to the Brunt- 
Vaisala frequency, N 

N' = --'/^^(^) (11) 
jp or or J 

which in our parameters looks like 

N2 = -/3p/3,^ (^Ij n2«-/52. (12) 

Here we explicitly wrote fip and fig to make clear that 
the Brunt- Vaisala frequency depends on the product of 
entropy and pressure gradient which can be different in 
global simulations. 
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1500 
time [27r/n„] 



2000 



2500 



Fig. 3. — In this run with = 1.0 a resolution of 576^ (run C, upper 
panel) and 1152^ (run D, lower palnel)and we turn off the entropy 
gradient after 800 local orbits (indicated by the black dashed line) and 
see how the instability decays. Enstrophy is shown with the black line 
and a-stresses with the blue line. Our fit is given through the red and 
green dashed lines respectively. We fit a decay time of t^z = — 1000 
for the enstrophy and = —400 for a. 



All quantities in Eq. (12) are positive. Thus the 
Brunt- Vaisala frequency is imaginary and therefore 
a linear buoyancy driven turbulence would have a 
amplification-rate T (x iN (x fi. However, we found that 
T oc fi^ provides a better fit. This once again reflects that 
the baroclinic vortex amplification is a non-linear effect. 
In linear convective instability a displaced parcel of gas 
feels a buoyancy force and thus accelerates propotion- 
ally to fi. But in the disk baroclinic instability first a vor- 
tex has to form with an azimuthal entropy gradient pro- 
portional to j6 (and T^ooi) and in a second step this vortex 
feels a torque proportional to j6. Therefore the amplifica- 
tion is proportional to ji^. The and \ooi dependance 
has also been derived by Lesur & Papaloizou (2010), see 
their Eq. (23) for an order of magnitude estimate of the 
growthrate. 

The amplification behavior in Fig. 1 already shows 
convergence for 576 grid cells resolution, e.g. 144/H in 
radial direction. 

If we compare our amplification timescales for the 
lowest entropy gradients with the migration times ob- 
tained by Paardekooper et al. (2010) we see that they 
are of the same order of magnitude. Which means that 
the vortex could have drifted into the central star be- 
fore it reaches strong a-values. However, Paardekooper 
et al. (2010) also state that their timescales refer to fully 
grown vortices of size H. Smaller vortices drift signif- 
icantly slower. This gives them enough time to reach 
a size, with which they provide sufficient angular mo- 



mentum transport, before they drift inward. 

To study the numerical dissipation effects even fur- 
ther we now assess how the vortices decay if baroclinic 
driving is switched off (Fig. 3). To do this we first evolve 
runs C and D with /3 = 1.0 and the two resolutions of 
576^ and 1152^ for 800 orbits and then turn off the en- 
tropy gradient so that /5 = 0.0. We observe that the vor- 
tices get smaller and that all relevant quantities like vor- 
ticity, co^, or a-stresses decay with exponential behav- 
ior. Godon & Livio (1999) saw the same exponential de- 
cay of vorticity when they analyzed longevity of anti- 
cyclonic vortices in protoplanetary disks. Their dissipa- 
tion was proportional to the effective viscosity applied 
in their numerical experiment. Here we find the same 
decay-rate for both resolutions, highlighting that the de- 
cay of vortices is no longer through numerical effects, 
but due to the radiation of waves as in Korotaev (1997). 

4.3. Saturation Values 

We have established that even shallow entropy gradi- 
ents lead to vortices but we still have to show that suffi- 
cient angular momentum transport can be reached with 
these shallow gradients. The saturation values of en- 
strophy, a>2, or Mrms are of interest as well. Note that 
we talk about saturation values of our 2D local simula- 
tions, where certain restrictions apply, see a more de- 
tailed discussion in the conclusions. In the next sec- 
tions we discuss the measured saturation values and an- 
alyze how the different controlling parameters influence 
amplification-phase and final values. 

4.3.1. Influence of Entropy Gradient 

In Fig. 4 we compare runs A, C and J (at a resolu- 
tion of 576^ and Tjiff = Tcqqi = 10) which differ only re- 
garding the value of /5. There is an initial exponential 
amplification-phase of a, E]^;n and lo^ that is shorter for 
high j6, followed by a saturated state. We also see that 
for lower /5 the saturation values are lower. We want 
to stress that we did not reach saturation for simula- 
tions J and K (at a resolution of 576^ and 1152^ and 
'''diff = "^cool = 10)- Even after 3000 local orbits vortex 
amplification was still ongoing. Here, Tjiff = 10 is much 
shorter than the amplification-rate we estimated in the 
previous section (t ~ 300). As we will see in the next 
section the amplification-phase is shortest if those time- 
scales are comparable, because T^j^ff also defines how 
fast pressure perturbations are damped. Although we 
expect the saturation values of simulation J and K to be 
higher than what they are right now, it is possible that 
they will still stay below the saturation values obtained 
in simulations with higher ^. 

The vorticity can be seen as a measure of the strength 
of the vortex. The higher the absolute value of the vor- 
ticity the stronger the vortex. The only stable vortices 

in disks are anticyclonic^ and therefore the vorticity has 
negative values. So the minimum value of vorticity 
(<^2,min) shows how strong a vortex is. To explain the be- 
havior of a;2,min (3rd panel in Fig. 4), cooling processes 
have to be taken into account. During the early phases 

^ Cyclonic vortices are also possible, but are quickly destroyed by 
shear (Godon & Livio 1999). 
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Fig. 4. — Time evolution of kinetic energy Ej.,,, (top), a-value (mid- 
dle) and minimum vorticity o^z^min (bottom) for a resolution of 576^ 
and Tjiff = Tjooi = 10 but different entropy gradientes: /S = 2.0 (green), 
^ = 1.0 (black) and ^ = 0.5 (red) (runs A, C, J). Saturation is first 
reached for high fi already after 300 orbits, then for /3 = 1.0 For /5 = 0.5 
no saturation is reached even after 3000 orbits. The increase in a;2,min 
after the point in time when saturation is reached can be explained 
through the heat transport across the vortex. Since it has reached its 
final and largest size heat transport takes longer due to the larger size 
of the vortex. 



thermalization is dominated by thermal diffusion (Pe- 
tersen et al. 2007b). As mentioned before this time-scale 
is shorter for smaller vortices. Therefore heat exchange 
between the vortex gas and the ambient gas is more ef- 
ficient than in later stages. Once the vortex has grown 
to its final size, thermal relaxation takes over. How- 
ever heat exchange in the center of the vortex is less effi- 
cient than in the earlier stages. The baroclinic feedback, 
e.g. the azimuthal entropy gradient across the vortex, is 
less efficient, the vortex grows weaker, and a;^ rises 
again, creating a flat yet extended vortex. 

4.3.2. Influence of Thermal Diffusion and Cooling Times 

We take a closer look at simulations with /5 = 1 and 
different combinations of K and Tcool to see how thermal 
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Fig. 5. — Comparison of different T^iff (right numbers) and Tcqqi (left 
numbers) for same /i = 1.0 (Runs C-I). The top panel shows the a-value 
and the bottom one !(rms. One can see that the early amplification- 
phase is determined by the diffusion time since the heating across the 
vortex is more important then vertical heat transport. We get faster 
amplification for higher T^^ff . Once the vortex grows larger heat trans- 
port gets more difficult and thermal relaxation dominates. Therefore 
the saturation values are determined through t^ooi. Saturation values 
are higher for shorter t^qoi 

diffusion and relaxation influence the saturation val- 
ues and the amplification-phases. As long as T^iff(;) = 

l^/K < Tcooi, T(jiff(/) will dominate the heat exchange 
from the inside of the vortex to the ambient disk. As 
the vortex grows Tjjff^;^ will increase and with that only 
contribute to the heat exchange at the outskirts of the 
vortex. T^ool will then dominate the interior of the vor- 
tex. 

For the simulations where we set Tjiff = TcqoI, TcqoI will 
take over when the vortex has reached a size of H. In 
radial extend this happens once the vortex has grown to 
its final size. 

This is consistent with what we see in Fig. 5. During 
the early amplification-phase simulations with equal 
Tdiff behave exactly the same. Eventually T^ooi takes over 
so that the saturation values are determined by t^oo\. For 
longer T^^oi saturation values are lower than for shorter 
^cool- 

4.3.3. Influence of Physical Domain 

A problem with local shearing sheet simulations is 
that eventually vortices grow to box-size. We cannot 
say whether they have reached their final size or just do 
not have any more room to grow. Another problem that 
arises with the periodic boundary conditions is that the 
vortices potentially interact with themselves and thus 
forcing (shaking) them to shed more waves and there- 
fore increase the a-values. To deal with that, we re-did 
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Fig. 6. — Snapshots of the z-component of the vorticity, Wz after 100, 
500, 1000, 1500 local orbits for the two different physical domains with 
/3 = 0.5. Initially both runs have vortices of equal size. Since there is 
less space between vortices, they can merge sooner in runs with the 
small physical domain. The vortices in the large physical domain take 
longer to grow. The dashed white box in the last plot indicates the 
area of the small physical domain. 



simulations A, C and J with a doubled physical domain 
(simulations A2, C2, J2 in Table 1). The resolution is 
the same. Instead of x = [—0.2,0.2] and y = [0.0,1.6] we 
switch to X = [-0.4,0.4] and y = [0.0,3.2]. We did not 
adjust the initial perturbation in any way. Therefore the 
initial state is perturbed at smaller wave numbers than 
in the smaller domain. If we go to even larger boxes the 
initial condition has to be adjusted so the the effective 
perturbation in the density is of the same strength as in 
the smaller physical domain. 

If we compare the time development of runs with a 
different physical domain (see Fig. 6), we see that vor- 
tices in fact do not merge as fast in the large domain be- 
cause there now is more space between them in radial 
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Fig. 7. — Time development of a and wl with fi = 0.5 for small (black) 
and large (red) physical domain (runs J and J2). Saturation values are 
lower in the large box than in the smaller box. 



direction, and they thus pass each other less frequently 
due to the extended azimuthal domain. Eventually they 
can merge as Godon & Livio (1999) saw, but the larger 
the box the longer it takes. We do not want to discuss the 
mechanism of how the process of vortex merging hap- 
pens exactly. This has been explained extensively in the 
field of fluid dynamics (see e.g. Cerretelli & Williamson 
2003). The merging process itself is not the focus of our 
study, because a) the vortex merging is strongly influ- 
enced by the box dimensions in a shearing sheet simu- 
lation and b) 2D flat vortices merge differently than full 
scale 3D vortices. The important thing is that vortices do 
indeed merge if the are sufficiently close to one another, 
but conserve co in the process. 

Another unphysical process that can occur in local pe- 
riodic simulations is that when the vortex approaches 
the integral scale it interacts with itself, the outer edges 
of the one side of the vortex almost touches the other 
side of the same vortex. We do not see this for the runs 
with the larger physical domain. Since the vortices in 
the larger domain do not interact with themselves, the 
saturation values are lower. However, they are still in 
the same order of magnitude (see Table 1). 

In Fig. 6 we show snapshots of the vorticity for /5 = 0.5 
(simulations J and J2). Initially there are several vortices. 
The larger ones sweep up the smaller vortices and thus 
grow further. At 150(5 local orbits there is only one vor- 
tex left for the small physical domain, whereas in the 
larger physical domain there are still three vortices. 

If we look at the a-value and enstrophy for these two 
simulations (see Fig. 7) we see that the value seems to 
decay in the larger box at the end of the run. How- 
ever this does not mean that the vortices die out. It 
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Fig. 8. — Vorticity profile (left) and a-stress (middle) for /S — 1.0 and the large physical domain (run C2). Yellow and red areas denote positive a- 
vaues whereas blue areas show negative a-stresses. In green areas a = 0. One can see the waves excited by the vortex. Those waves are responsible 
for the angular momentum transport. It is a localized process. Since the vortex and the vorticity-waves fill out a smaller area of the box in the 
large box (large green areas where there is no angular momentum transport) and our calculation of the saturation values averages over the entire 
area of the box, the saturation values seem to be lower The plot in the right upper panel shows an azimuthal average over the UxUy. Inside an 
ideal vortex a-stresses would sum up to zero. However, as indicated in the lower right plot, the vortex has a complex structure which leads to 
deviations from the idealized shape. 



more so reflects fluctuations in the vortex interaction, 
modulating a, as also can be seen in the small domain 
case at high frequency. We calculate the values as a 
mean over the entire box but especially the angular mo- 
mentum transport is a very localized process as can be 
seen in Fig. 8 (this time for /5 = 1.0 after 1000 orbits). 
Here we show the product UxUy at each location in the 
box. Most areas of the box have an M^Wy-value close 
to zero. However, one can clearly see bands excited 
by the vortex with positive UxWy-values. These bands 
are inertia-acoustic waves which are responsible for the 
angular momentum transport (Klahr & Bodenheimer 
2003; Mamatsashvili & Chagelishvili 2007; Heinemann 
& Papaloizou 2009; Tevzadze et al. 2010). If we had an 
ideal vortex with a smooth surface we would expect that 
Ml- My sums up to zero within the vortex. However the 
vortex has a more complex structure as can be seen in 
the lower right plot of Fig. 8. This leads to an negative 



net a-value across the vortex. 

To properly compare the values of a for both physical 
domains, the box average has to be taken. If the average 
over an equal physical size centered around a vortex, as 
indicated by the white dashed lines in Fig. 6, is taken, 
then the a-values agree again. The a-values are gener- 
ated only in the vicinity of vortices. 

4.4. Correlations 
It is a feature of baroclinic instability that the satu- 
ration values of Urms/ Prms Seem to correlate with 
each other. In Fig. 9 we plot the different quantities as 
a function of ol. Figure 9 shows the dependencies on a 
for all our simulations. The colors represent the differ- 
ent entropy gradients: /5 = 2.0 (black), /5 = 1.0 (red) and 
/5 = 0.5 (green). The different combinations of diffusion 
and cooling times are represented through the different 
symbols. We find that the following relations are good 
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Fig. 9. — Saturation values of <jj\, prms and linns as a function of sat- 
urated (value at the end of the simulation for f> = 0.5) a-value and 
all our runs with the small physical domain (runs A-P). The symbols 
show the different combinations of t^oqi (left numbers) and T;ijff (right 
numbers). Where red are runs with black /5 = 2.0, p = 1.0 and green 
/5 = 0.5. The black dashed line shows the dependency that we fit. 

fits to our simulation results 

Mrms=3yaCs (13) 
Prms =2^^,00 (14) 

<x;2 = 5ang. (15) 

We can derive the typical length-scale of angular mo- 
mentum transport L, of the system if Eq. (13) is inserted 
into the general a formalisms (Shakura & Sunyaev 1973) 

V = CiCgH = Mrmsi' SO that 



(16) 



indicating smaller structures than the vortices in our 
simulations and also smaller than the vorticity in stan- 
dard a-models where cv oc -^a with a different coefficient 
(Cuzzi et al. 1994). 
We do not perform a more exact analysis of these 



10" 



10" 



X 



10-^ ^ - - 



10" 



+ x,^f10 
A x,^,=30 
X x,^,=100 

• x„f10 

• x,„f30 

• x„f100 



0.1 



1.0 



10.0 
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depending on /S. Runs with parentheses around them were not satu- 
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account when we fit the a — /i-relation.The symbols show the different 
combinations of t^qoI (symbols) and T^iff (colors). 

dependencies (varying initial conditions) before we do 
three-dimensional simulations. 

4.5. Dependence on /? 

In Section 4.2 we showed that amplification of vortices 
for low entropy gradients is computationally demand- 
ing in terms of evolution time. Thus it is difficult to ex- 
tract saturation values for entropy gradients even shal- 
lower than [i — 0.5 with the computational resources at 
hand. 

In Fig. 10 we plot the a-stresses as a function of the 
entropy gradient. Note that we choose a different color- 
coding than in Fig. 9. Here symbols represent the ther- 
mal cooling times whereas colors represent thermal dif- 
fusion times. The dashed black line illustrates a slope 
CK j]P-^ which is a reasonable fit for the set of points 
with Tpooi = 30,Tjjjff = 10 (black triangles) and = 
100,Tdiff = 30 (orange x). We cannot predict a-values for 
specific entropy gradients and thermal cooling and re- 
laxation times. 

The key issue is less a strong correlation between a 
and f> but rather the lack thereof. The strength of the a- 
stresses reflects the size and the amplitude of the largest 
vortex. Its size is defined by H only and not by any 
of the other t and /5 parameters. As long as t and 
/S are sufficient to replenish vorticity at the loss-rate, 
the a-stresses should be independent of t and j6. The 
loss time-scale via generation of waves and Re3molds 
stresses is rather long, see Section 4.2 and Fig. 3. Thus 
as long as the amplification-rates are faster than decay- 
rates one should always obtain roughly the same a- 
values. 

5. SUMMARY AND CONCLUSION 

In this paper we have conducted an extensive param- 
eter analysis for the baroclinic vortex amplification. In 
particular we analyzed the influence of the global en- 
tropy gradient, thermal relaxation and cooling as well as 
numerical parameters such as resolution, box size, and 
amplification-rates for vortices and saturation values of 
a. 

The most important result of our study is that we find 
vortex growth even for entropy gradients as low as /3 = 
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0.5. However the amplification rate is of the order of 
several 100 local orbits which makes it difficult to extract 
reliable saturation values for the efficiency of angular 
momentum transport. 

Recently Paardekooper et al. (2010) studied the migra- 
tion behavior of vortices in global accretion disks. They 
found significant radial drift for fully grown vortices 
with drift times shorter than the vortex amplification 
times we measure in this paper. Nevertheless, this is not 
a contradiction, because as also shown in Paardekooper 
et al. (2010) drift rates strongly depend on vortex size. 
Thus the tjrpical life cycle of a growing vortex might 
be starting as a growing small vortex without relevant 
radial drift, which starts drifting as soon as it reaches 
its saturated state. Therefore radial drift does not affect 
the study of vortex amplification discussed here. How- 
ever, it will affect the time a single vortex can partake in 
angular momentum transport. Future work will have 
to investigate radial drift of growing vortices in global 
simulations. Note here that Paardekooper et al. (2010) 
studied the migration in barotropic disks, in which no 
vortex amplification occurs. 

The amplification-phase of the vorticies can be mea- 
sured in the strength of the overall velocity fluctuation 
which seem to be growing exponentially on a certain 
time-scale t oc |S~^. Therefore amplification for steeper 
entropy gradients is faster, i.e. t = 16 for j6 = 2.0 and 
T — 70 for |S — 1.0. With these short amplification-times 
we do reach saturation. Whereas the ^ = 0.5 was still 
growing after 3000 orbital periods, when we stopped 
tihe simulation. 

Other parameters that influence the evolution of a- 
stresses are the thermal cooling and relaxation times. 
The diffusion times define the amplification phase of the 
vortices because diffusion dominates small scales, e.g. 
small vortices. We see faster amplification for longer 
diffusion times. Cooling time on the other hand deter- 
mines the saturation values. Here, longer time-scales 
produce lower saturation values. 

For the angular momentum transport we get a-values 
up to 10^2 for ^ ^ 2.0 and lO"^ for ,6 = 1.0 and ^ = 0.5. 
These values are not so different to the ones found with 
MRI in active layers (Flock et al. 2011) and stronger than 
the 10~^ found in dead zones (Dzyurkevich et al. 2010), 
which shows that entropy gradients can be an impor- 
tant mechanism to transport angular momentum in a 
dead-zone. Realistic entropy gradients in protoplane- 
tary disks are around /S — 0.5 and ^ — 1.0 which can 
be derived out of the data obtained by Andrews et al. 
(2009) as discussed in Klahr (2013 submitted to ApJ). Al- 
though we could not reach saturation in all our simula- 
tions for these entropy gradients we do see reasonable 
a-stresses of the order of 10~^ to 10"^. We expect the 
final values to be in this range which still provides suf- 
ficient angular momentum transport in a disk. Yet, we 
have to consider certain cavities: 1.) Our simulations are 
2D simulations and lack the 3 dimensional structure of 
the vortices. This might very well affect the strength of 
the a-values. 2.) We do not consider migration of vor- 
tices, but rather have periodic boundary conditions. It 
is not clear for how long vortices can play a role in an- 
gular momentum transport before they migrate into the 



central star. Thus we cannot say how many vortices are 
in a disk at any given time. The higher the number of 
vortices, the higher the a-values will be. The interplay 
between migration and Rejmolds stresses definitely has 
to be analyzed in future models. 3.) The formation pro- 
cess for vortices is still not clear. It is unknown how long 
the initial formation of a vortex takes, by which process 
they are formed and if there are processes which can de- 
stroy them before the reach full growth. Therefore, our 
saturation values have to be viewed with caution and 
cannot be seen as face values for protoplanetary accre- 
tion disks. As relation between entropy gradient and 
strength of angular momentum transport we only find 

a weak dependence of a (x ji^^^. 

Since local simulations are always limited by the box 
size we also conduct simulations in larger boxes. We do 
not see a difference in the initial amplification-phase. At 
later stages the amplification last longer for larger boxes 
and also is slower. Since part of the vortex evolution 
happens through merging of smaller vortices, growth 
takes longer in larger boxes simply because there the ra- 
dial distance between vortices is bigger and thus merg- 
ers are less likely. 

The saturation values of velocity fluctuations reached 
for the larger box sizes are slightly lower than for the 
smaller box sizes. This is due to two reasons. One is that 
we see some artificial enhancement in vortex strength 
in the smaller box. Once the vortex has reached box- 
size it can no longer grow. It is forced to interact with 
itself thus emitting more waves. This does not happen 
in larger boxes. 

The other reason is that the number of vortices per ra- 
dial distance is independent of box size because their 
typical maximum size is in the order of a pressure scale- 
height. In the azimuthal direction the number of vor- 
tices is limited to 1 per radius, because otherwise merg- 
ing will occur on short time-scales. Therefore the overall 
density of vortices per simulation volume (area) is lower 
in simulations with the larger azimuthal extend. Here 
we want to note that our larger boxes with H/ r — 0.1 
and Ly — 32 are only a factor of about two shy of the 
equivalent In global simulation. 

Overall, we conclude that the baroclinic vortex am- 
plification works reasonably well for entropy gradients 
as low as /3 = 0.5. This /3 corresponds to a Richardson- 
number of Ri = —1.5 X 10^'^. This makes BVA a rel- 
evant mechanism for angular momentvim transport in 
the dead-zone. 

An exploration of lower entropy values will have to 
be postponed due to the long evolution time required. 
In the future we will study stratified 3D boxes and the 
interaction of dust with the vortices. 

Our simulations were conducted partly on the MPIA 
cluster THEO in Garching, and on the JUGENE machine 
of the JSC using the grand HHD19. This work was par- 
tially supported by the National Institute for Computa- 
tional Sciences (NICS) under TG-MCA99S024 and uti- 
lized the NICS Kraken system. This collaboration was 
made possible through the support of the Annette Kade 
Graduate Student Fellowship Program at the American 
Museirai of Natural History. NR also wants to thank 
IMPRS-HD. 
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the identical behavior for the isothermal case (dashed and dotted lines), yet in the adiabatic case the TRAMP code shows an artificial amplification 
of kinetic energy (dashed-dotted line). The Pencil Code does not show this behavior. 



APPENDIX 
NUMERICAL ARTEFACTS 

Shearing sheet simulations with the TRAMP code have displayed unreliable behavior for the extreme cases of 
cooling times, either isothermal (TcqoI — 0) or adiabatic (TcqoI — °°)- In the first case, a global pressure gradient in a 
locally isothermal disk leads to the amplification of radially propagating soimd waves, which is a physically realistic 
case (see the derivation in Klahr 2013 ApJ submitted), but only shows up in local radially periodic simulations because 
the sound wave can propagate through the the box for an unlimited amount of time, which of course is not possible 
in a global disk. This physical instability can thus be found both in ID radial TRAMP as well as in PENCIL CODE 
simiilations with remarkably identical growth behavior. This means, having a too short cooling time artifacts from 
these radially propagating sound waves coiild ruin our models. Nevertheless, as pointed out by Klahr (2013 ApJ 
submitted) already a cooling time of Tcooi — 0.01 will suppress these sound wave instability completely. 

On the other hand the adiabatic simulations using the TRAMP code were showing a weak amplification of kinetic 
energy over very long time scales which is the accumulation of numerical error in the quasi dissipation free TRAMP 
scheme. This behavior is independent of the chosen entropy gradient and results from the conservative treatment 
of Coriolis forces. Again the PENCIL CODE with its explicit dissipation does not allow for this accumulation of this 
numerical error, even in the presence of a radial entropy gradient (see solid and dashed-dotted line in Fig. 11). 
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